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Abstract 

Background: Many mathematical models characterizing mechanisms of cell fate decisions have been constructed 
recently. Their further study may be impossible without development of methods of model composition, which is 
complicated by the fact that several models describing the same processes could use different reaction chains or 
incomparable sets of parameters. Detailed models not supported by sufficient volume of experimental data suffer 
from non-unique choice of parameter values, non-reproducible results, and difficulty of analysis. Thus, it is necessary 
to reduce existing models to identify key elements determining their dynamics, and it is also required to design the 
methods allowing us to combine them. 

Results: Here we propose a new approach to model composition, based on reducing several models to the same 
level of complexity and subsequent combining them together. Firstly, we suggest a set of model reduction tools 
that can be systematically applied to a given model. Secondly, we suggest a notion of a minimal complexity 
model. This model is the simplest one that can be obtained from the original model using these tools and still able 
to approximate experimental data. Thirdly, we propose a strategy for composing the reduced models together. 
Connection with the detailed model is preserved, which can be advantageous in some applications. A toolbox for 
model reduction and composition has been implemented as part of the BioUML software and tested on the 
example of integrating two previously published models of the CD95 (APO-1/Fas) signaling pathways. We show 
that the reduced models lead to the same dynamical behavior of observable species and the same predictions as 
in the precursor models. The composite model is able to recapitulate several experimental datasets which were 
used by the authors of the original models to calibrate them separately, but also has new dynamical properties. 

Conclusion: Model complexity should be comparable to the complexity of the data used to train the model. 
Systematic application of model reduction methods allows implementing this modeling principle and finding 
models of minimal complexity compatible with the data. Combining such models is much easier than of precursor 
models and leads to new model properties and predictions. 



Background 

Systems biology aims to study complex interactions in 
living systems and focuses on analysis and modeling 
their properties. Mathematical modeling provides several 
ways to describe biological processes based on experi- 
mental information of different kind. However, creation 
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of detailed models not supported by enough experimen- 
tal data often makes their analysis and interpretation dif- 
ficult [1]. Several aspects of the same process can be 
modeled using different levels of abstraction involving 
different reaction chains, chemical kinetics, and incom- 
parable sets of parameters. Such models are difficult to 
merge. Meanwhile, merging is an important approach 
for creation of complex models. Thus, development of 
efficient methods and software, allowing us to combine 
models, is the object of intense study in systems biology. 
In our work, we focused on the fact that, generally, com- 
plexity of models is not comparable to the volume of 
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experimental data used to adjust their parameters. Due 
to this fact, we turn to the methods of model reduction 
allowing us to minimize models complexity without af- 
fecting the model simulation dynamics. 

Model reduction is a well-established technique in 
many fields of biochemical research and engineering. It 
has been used for many years in chemical kinetics (for 
reviews, see [2-4]) and has already found multiple appli- 
cations in systems biology, including discrete modeling 
[5] and modeling of metabolic pathways [6,7]. The prin- 
ciples of this technique have been applied in computa- 
tional biology [8] and implemented as a part of widely 
used pathway simulators such as BioUML [9], BIOCHAM 
[10], COPASI [11], and GINsim [12]. Model reduction led 
to new insights in mechanisms of translation regulation 
by microRNAs [13,14] and was applied for analysis of 
such signaling pathways as JAK-STAT [15], NF-kB [16], 
and EGFR [17]. 

In our investigation, we used the principles of model 
reduction to construct reasonably accurate minimal size 
approximations of two different models describing the 
CD95 signaling pathways [18,19]. The first model ex- 
plores pro-apoptotic properties of CD95 after stimu- 
lation by its natural ligand CD95L or by agonistic 
antibodies anti-CD95 implying formation of the death- 
inducing signaling complex (DISC) [18,20]. DISC con- 
sists of oligomerized CD95, death domain-containing 
adaptor FAS-associated molecule (FADD), procaspases-8 
and -10, and two isoforms of cellular FLIP (CFLAR) 
protein (cFLIP long and cFLIP short). Caspase-8 leads to 
activation of effector caspase-3 directly (type I cells) or 
via stimulation of cytochrome C release from the mito- 
chondria (type II cells) [21]. The latter step requires for- 
mation of the apoptosome complex and activation of 
caspase-9. Once activated, caspase-3 cleaves poly(ADP- 
ribose) polymerase (PARP), thereby making the apop- 
totic process irreversible. The second model describes 
the state when CD95 not only activates the pro-apoptotic 
pathway, but also induces transcription factor NF-kB that 
is an important regulator of cell survival functioning [19]. 
This is possible due to cFLIPL cleavage in the DISC com- 
plex. The cleaved p43-FLIP directly interacts with the IKK 
complex and activates it. The activated IKK performs 
phosphorylation of the IkB protein and thereby frees 
NF-kB. 

The authors of the models have evaluated concen- 
tration changes of major apoptotic molecules by 
Western blot analysis and represented them as rela- 
tive values at given time points. Using the systematic 
methodology [2] implemented in the BioUML soft- 
ware, we reduced the models so that they still satisfy 
these data. This allowed us to simplify the overlap- 
ping components of the models and find a way for 
their composition. 



Methods 

Model reduction 

Mathematical modeling of biological processes based on 
the classical theory of chemical kinetics assumes that a 
model consists of a set of species S = (S lt ._S m ) associ- 
ated with a set of variables C(t) = (Ci(£), . . ., C m (t)) de- 
pending on time te [0,T\,Te R + and representing their 
concentrations, and a set of biochemical reactions with 
rates v(t) = (vi(t), . . ., v n {t)) depending on a set of kinetic 
constants K. Reaction rates are modeled by mass-action 
or Michaelis-Menten kinetics. A system of ordinary dif- 
ferential equations (ODE) representing a linear combin- 
ation of reaction rates is used to describe the model 
behavior over time: 

^l = N-v(C,K,t),C(0) = C°. (1) 

Here N is a stoichiometric matrix of n by m. We say 
that C ss is a steady state of the system (1) if 

N-v(C ss ,K, t) = 0, lim Q(t) = Cf . (2) 

t^oo 

Model reduction implies transformation of the ODE 
system to another one with smaller number of equations 
without affecting dynamics of variables Ci(t), . . ., C s (t), 
which is fixed by a set of experimental points Cf^fey) at 
given times t^j = 1, . . ., r if where r t is number of such 
points for the concentration Q(£), i = 1, . . To check 
dynamics preservation in the course of model reduction, 
we consider the function of deviations defined as a nor- 
malized sum of squared differences [11]: 

/-*(<?,*) = ±±°^{Q( tlj ) - CT(%)) 2 , 

i=l /=1 1 

(3) 

where o min = min 0)i, weights Oi = yj rj l • ^ [Cf p (%) ) 2 

are calculated as mean squared values of experimental 
concentrations for each variable and the normalizing fac- 
tor (x) min /(x>i is used to make all concentration trajectories 
to have similar importance. 

Reducing kinetic model is possible when some quan- 
tities are much smaller than other quantities and can be 
neglected. Usually, this implies some qualitative relations 
(much bigger, much smaller) between model parameters. 
When these relations satisfy certain rules, we can ap- 
proximate the detailed model by a simpler one. If the 
parameters are (approximately) determined, as in our 
case, we find an approximation specific for a region of 
the parameter space. However, if we want to investigate 
the model behavior for the entire space, we need to de- 
compose it into regions characterized by asymptotically 
different behaviors of the dynamical systems. Afterwards, 
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a specific reduction should be performed for each region 
(Figure 1). 

We define the minimal complexity model as a reduced 
model with the minimal number of elements (species 
and reactions) so that the deviation function value (3) 
calculated for the reduced model is different from the 
original model no more than 20%. Such threshold is 
explained by the fact that in this work we considered ex- 
perimental data obtained by Bentele et al. [18] and Neu- 
mann et al. [18] using Western blot technology with the 
standard deviation 15-20%. 

Reduction of mathematical model complexity is 
achievable by different methods [2]. Description of the 
methods directly used in our work is provided below. 

(MR1) Removal of slow reactions. We say that reaction 
ri is much slower than reaction r 2 , v n (t)«v r2 (t) if 
max \v n (t)\ < k - max \v r2 (t)\ , where k=10' 2 . When such 

reactions do not affect experimental dynamics of species, 
we neglect them. Note, that in some cases, it is sufficient 
to consider k- 10 _1 or even k- 1. 

(MR2) Quasi-steady -state approximation [22] assumes 
that all variables C(t) of the model are split in two 
groups: basic variables C s (t) and quasi-stationary vari- 
ables Cf [t) for which we could write: 



dC s (t) 
dt 



where e is small parameter. We consider the approxima- 
tion £— >0 resulting in the equation • v^(C, K, t) = 0. 
Thus, the chain of reactions S— >/— >P with a quasi- 
stationary species / can be reduced to the form S — > P. 



(MR3) Lumping analysis refers to reducing the num- 
ber of model variables by grouping some species. In par- 
ticular, this can be illustrated by reactions 



rl : A + S^P^ r2:A + S 2 



(4) 



with kinetic rates v r i(t) — /ti-Ca(^)-C < s 1 (^) and v r2 (t) = 
k 2 <C A (t).C S2 (t). lik 1 = k 2 and C Sl (t) = C S2 {t), then the 
system (4) can be replaced by a single reaction 

2A + 

where C s (t) = C Sl (t) = Cs 2 (t)> C P {t) = C Pl (t) = Cp 2 (t) 
and k = ki = k 2 . Note, that in this example lumped vari- 
ables Cs x (i) , Cs 2 (t) and C Pl (t) , Cp 2 (t) are linearly 
dependent: however, this is not required in general case. 
For conditions on lumpability in monomolecular reac- 
tion networks, see [23]. 

(MR4) Removal of approximately linearly dependent 
variables. If variables A(t) and B(t) (species concentra- 
tions or reaction rates of the model) are approximately 
linearly dependent: 

A(t)*k-B(t), te[0,T\, k = const, 

then we can replace one of them with another in all kin- 
etic laws of the model. 

(MR5) Simplification of the Michaelis-Menten kinetics. 
Consider a reaction r of the form S -E — > P, where 
an enzyme E converts a substrate S into product P. 



C A (x) »C B (t) 



[ 



C A (t) » C,(t) v rl (t) ^ v r2 (.t) 



4 



Complete mode! 



Parameter space] 



| - w 

V rl 0) * y r2 f» 



C ] H3>»[~P 



CM » cAt) 



V-LH*4 D 

-I— 



T 



r ^ 



C B (.t) » C A (t) v rl (t) w v r2 (t) 



tf rl (0 * v r2 (0 



Figure 1 Schematic representation of the model reduction technique. A model consisting of four species {A, B, C, and D) and two reactions 
(r1 and rl) was reduced taking into account the relations between time-dependent concentrations C A {t) and C fi (f), as well as reaction rates v n {f) 
and v r2 {t). In the case when C A {t) » C B {t), the concentration of A can be considered constant for some period of time. The same simplification 
could be applied to B. If v n {t) ~ v r2 {t), then C could be removed, resulting in direct formation of D triggered by the binding of A and B only. 
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Reaction rate of such reaction is frequently defined by 
the Michaelis-Menten formula 



V r {t) 



k.c s (tyc E (t) 

Km + C s {t) 



(5) 



where the enzyme concentration is a dynamic variable 
allowing to use the same kinetics in different regions 
of the phase space. If Km » C s (t), then we can reduce 



this formula to the form: v r (t) 



k.C s (t)-C E (t) 
Km 



. On the other 



hand, if C s {t) » Km, then (5) can be replaced by the 
equation v r {t) = k • C E (t). 

(MR6) Simplification of equations based on the law of 
mass action when one reactant dominates others. Con- 
sider a reaction r of the form Si + S 2 — > P with the kin- 
etic law 



v r (t)=k-C Sl (t)-C S2 (t)- 



(6) 



When C Sl (t)»C S2 (t) for ts [0,7], the formula (6) can 
be replaced by the linear equation v r (t) = 
k-Cs 1 (0)-Cs 2 (t) . The time T of validity of such pseudo- 
lineary approximation is defined as a period during 
which the relative change of Cg^t) does not exceed 

L.e. k I Cs 2 (t)dt < e. 



'I 

J o 



One way to perform the model reduction is to apply 
the foregoing methods in the numerical order (Figure 2). 
Note that the methods MR5 and MR6 are of the same 
type, so it does not matter which one to use first. This 
way assumes that we calculate the value of the distance 
function (3) at each reduction step to test whether ap- 
proximation of the experimental data is still within al- 
lowable limit or not. 



I Removal of slow reactions I 



Qua si -steady-state 
approximation 



3 

I Lumping analysis I 

r 



Removal of linearly 
dependent variables 



T 



Simplification of 
equations with the 
Micha elts-Menten kinetics 



Simplification of 
equations based on the 
law of mass action 



Figure 2 Flow chart of the model reduction. The methods of 
model reduction are presented in the order of their application. The 
last two methods are of the same type, so this is a choice of the 
systems biologist which one to use first. 



Model analysis and comparison 

Model comparison using Akaike information criterion 
(AIC). This criterion [24] defines the relative complexity 
of a model based on the goodness of experimental data 
fit and the number of model parameters \K\ (initial con- 
centrations of species are not considered): 



AIC=x* + 2-\K\. 



(?) 



The function x is defined by the formula (3) with 
weights l/o| (instead of cd m i n lo)j) [15], where mean devi- 
ations Oij of experimental values Cf cp (tij) are calculated 
using smoothing spline [25]: 



4=iE(^(^)- c r(^)) 2 - 

*k=-2 



For ; + k < 0 and / + k> r t we assumed Q(^y +y ^) = Cf*^ 

(k J+ k) = o. 

Models can be compared using the Akaike criterion, if 
they approximate the same experimental data. When the 
AIC value of one model is less than of others, we say 
that this model is simpler in terms of this criterion. 

Model comparison with the mean AIC. When we want 
to compare models approximating different sets of ex- 
perimental data, we could calculate the mean AIC coeffi- 
cients by the formula 



AIC m 



AIC 



l exp 



(8) 



where n exp determines the number of experimental 
points. 

Steady-state analysis finds values of species concentra- 
tions according to the rules (2). It is desirable to reduce 
a model so that steady-state concentrations of all experi- 
mentally measured chemical species have not changed 
significantly. 

Sensitivity analysis reveals steady-state concentrations 
response to parameter perturbations. The local sensitiv- 
ities Cf , /=1,..., m, for all parameters p^ j = 1, . . ., \K\, 
are calculated via finite difference approximations: 

c _8Cf _C?(pj+Apj)-C?(pj) 



dpj 



A Pj 



It is useful to compute the mean log sensitivity for all 
fitted species: 



£ iV ln N 



%*0, i= j=l,...,n p , (9) 



where n p is the number of all parameters \K\ or the num- 
ber of parameters retained after the model reduction. 
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Modular modeling 

Considering the mathematical models of CD95 signaling 
pathways [18,19], we decomposed them into functional 
modules. This step allowed us to identify overlapping 
components of the models and simplify their analysis. 
We defined a module as a submodel (including several 
species, reactions and parameters of the model) with in- 
put, output and contact ports. The first two types of 
ports characterized variables calculated in one module 
and passed to another through a directed connection. 
The contact ports declared common variables of mod- 
ules via undirected connections (for more details of 
modular modeling, see [26]). 

Results 

Preliminary analysis, problems and inconsistencies in the 
precursor models 

The mathematical model by Bentele et al (Figure 3A) con- 
sists of 43 species including a virtual variable x aa intended 
for quantification of "apoptotic activity" caused by active 
caspases-3, -6 and -7 and calculated by the formula 

dx aa (1—Xaa)' {^36act ' Qasp3 + ^36act'C C asp6 ~l~ ^7 'act *C ' caspl) 7 

dt Km 367act + (1 

Xaa ) 

(10) 

where k 36act , k 7act , Km 367act e K. The variable x aa specifies 
a process of degradation in the model introduced as an ex- 
ponential decay function fdegX^aa)- This function is defined 
by the formula • x\ a + 

kds with k d Mds G K for all active 
caspases and complexes containing their cleaved products, 
and by the formula ka • x aa for all other molecules besides 
cPARP (cleaved PARP), where fdegX^aa) is constant. All 
species in the model are degraded with the exception of 
cytochrome C and Smac stored in the mitochondria, 
whose concentrations are constant. 

In total, the model contains 80 reactions (Table 1) includ- 
ing 24 reactions based on mass action kinetics, 12 reactions 
taken with kinetics of Michaelis-Menten, 41 reactions of 
degradation, the reaction of x aa production specified above 
and two reactions of cytochrome C and Smac release from 
the mitochondria modeled using a discrete event. The lat- 
ter implies complete cytochrome C and Smac release 
within 7 minutes as soon as tBid reaches a certain level in 
comparison to Bcl-2/Bcl-XL. Since the authors did not pro- 
vide exact form of the release function, we proposed the 
sigmoid function which has the expected behavior: 



1 H~ €Xp(j< con t r ' ( t + t trigger H~ 0 '-5 t release)^) 

(ii) 

where t trigger is the time when tBid concentration reaches a 
value of Bcl-2/Bcl-XL, t re i ease is the start time of release and 
Kontr is the contraction coefficient. 



The model comprises 43 species (including x aa ) and 
45 kinetic parameters, which the authors estimated 
based on the experimental data obtained by Western 
blot analysis for the human cell line SKW 6.4. Cells were 
stimulated by 5 ug/ml and 200 ng/ml of anti-CD95 
(fast and reduced activation scenarios, respectively) 
and dynamics of several proteins (Bid, tBid, PARP, 
cPARP, procaspases-2, -3, -7, -8, -9, cleaved pro- 
duct of procaspase-8 p43/p41, and caspases-8) were 
measured. 

We could not reproduce the dynamics of the original 
model using the parameter values provided by the au- 
thors. In particular, there was too rapid consumption of 
procaspases-2, -3, -7, -8 in the case of 5 ug/ml of anti- 
CD95 and procaspases-2, -9 in the case of 200 ng/ml. 
Degradation rates of procaspase-8 and caspases-8 was 
insufficient for both activation scenarios. Thus, we had 
to make several modifications of the original model to 
obtain the same dynamics as described in the original 
paper. Namely, we multiplied the rate constants of all 
the bimolecular reactions by the value 5- 10" 5 and speci- 
fied k d equal to 3.56 min 1 and 0.62 min 1 (instead of 
0.891 min" 1 and 0.184 min 1 ) for fast and reduced activa- 
tion scenarios, respectively. 

The model by Neumann et al includes 23 species, 23 
reactions constructed according to the law of mass ac- 
tion, and 17 kinetic parameters (Figure 4A). It repro- 
duces experimental measurements of 8 entities (total 
amount of procaspase-8, its cleaved product p43/p41, 
caspase-8, procaspase-3, caspase-3, Ii<B-a, phosphory- 
lated Ii<B-a and cleaved p43-FLIP) obtained using West- 
ern blot technology for HeLa cells stably overexpressing 
CD95-GFP and treated with three different concentra- 
tions of agonistic anti-CD95 antibodies (1500, 500 and 
250 ng/ml). This is the reduced model suggested by the 
authors for description of the CD95 -mediated pathways 
of cell death and NF-kB activation. Model reduction was 
based on intuition about what is important in the bio- 
chemical mechanisms. It was performed as long as the 
parameter estimation procedure (repeated after each 
step of the reduction) was able to provide a good fit of 
the experimental data. Particularly, the details of some 
complex formations were omitted as well as some non- 
measured species were removed. The authors submitted 
the model to the BioModels database [27]. Thus, we had 
no problems reproducing it. 

Reduction of the CD95-signaling model 

We started the Bentele s model reduction by excluding 
the direct dependence of the virtual species x aa on 
caspases-6 and -7. For this purpose, we approximated 
the amount of caspases-7 by the linear function of 
caspase-3 concentration: 
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Figure 3 The model by Bentele et al. and results of its reduction. A. The original model decomposed into modules according to three steps 
of apoptosis: activation of caspases-8 and -9 and inactivation of PARP. The species retained after the model reduction are colored. B. The modular 
view of the model. The dashed connections were deleted during the model reduction. C. The minimal reduced model. D. The graphical notation 
used for representation of the models A-C. 
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Table 1 Summary of reactions from the original model by Bentele et al. and the reduced model 



brl CD95L + CD95R -> CD95R:CD95L (MA) 

br2 FADD + CD95R:CD95L -> DISC (MA) 

br3 pro8 + DISC -> DISC:pro8 (MA) 

br4 pro8 + DISC:pro8 -> DISC:pro8 2 (MA) 

br5 DISC:pro8 2 -> DISC:p43/p41 (MA) 

br6 DISC:p43/p41 -> casp8 + DISC (MA) 
Activation of caspase-8 by caspase-3 



Original model 


Reduced model 


Me Reactions (Kinetics) 


Rates Reactions (Kinetics) 




(nM/min) 


Activation of caspase-8 induced by CD95 



071 0 U brl * CD95L DISC (k LR ■ C CD95R ■ C CD95Ll C CD95R 

oVi o° is ^ eterminec ^ by tne fo rmu ' a 04)) 



0°/1 0" 1 br2- pro8 -DISC - casp8 (k Dlsc ^ ro8 ■ C pro8 ■ 
0710" 1 Cdisc) 



br7 pro6 -casp3 — ► casp6 (M-M) 
br8 pro8 -casp6 -> casp8 (M-M) 
Inhibition of the DISC complex 



071 0" 1 br3* pro8 -casp3 casp8 (k 38 ■ C casp3 ) 



br9 DISC + cFLIPL -> DISC:cFLIPL (MA) 

br10 pro8 + DISC:cFLIPL -> DISC:cFLIPL:pro8 (MA) 

brl 1 DISC + cFLIPS -> blocked DISC (MA) 

br12 DISC:pro8 + cFLIPL -> DISC:cFLIPL:pro8 (MA) 

br13 DISC:pro8 + cFLIPS -> blocked DISC (MA) 

br14 DISC:cFLIPL:pro8 -> p43/p41 + blocked DISC (MA) 

Activation of caspase-9 triggered by cytochrome C 



0710° br4* cFLIP + 2- DISC + pro8 -> DISC:FLIP:pro8 

n o (k D isc_FUP ■ Cflip ■ Cdisc) 



0" 2 /10" 3 br5* DISC:FLIP:pro8 -> p43/p41 



(^DFp8 ■ C D isc.FLIP:pro8) 



brl 5 Cyt C Stored -> Cyt C {f release® -CcytCstored) 

br16 Apaf-1 + Cyt C -> Cyt C:Apaf-1 (MA) 
brl 7 Cyt C:Apaf-1 -> Apaf-1 + Cyt C (MA) 
brl 8 pro9 + Cyt OApaf -> Apop (MA) 
brl 9 Apop -casp3 -> casp9 + Cyt C:Apaf-1 
(M-M) 

br20 Apop -> casp9 + Cyt C:Apaf-1 (MA) 
Activation of Bid 



br6* pro9 — > cleavage (the formula (16)) 



br21 pro2 -casp3 — > casp2 (M-M) 
br22 Bid -casp8 tBid (M-M) 
br23 Bid -casp2 -> tBid (M-M) 
Blocking of IAP by Smac 



br7* Bid -casp8 - tBid (^ ■ C casp8 ■ C Bid ) 



0"71 0° br8* pro2 -casp3 -> cleavage (^g^) 



br24 Smac stored -> Smac (/ > e/eose (f) '-C Smac stored ) 
br25 Smac + IAP -> IAP:Smac (MA) 
br26 IAP:Smac -> Smac + IAP (MA) 
Activation of caspases-3 and -7 



br27 pro3 -casp8 casp3 (M-M) 
br28 pro3 -casp9 casp3 (M-M) 
br29 pro7 -casp8 -> casp7 (M-M) 
br30 pro7 -casp9 — > casp7 (M-M) 
inactivation 



0 71 0 U br9* pro3 -casp8 casp3 (M-M) 
0" 3 /10° 

0°/10" 1 br10* pro7-casp8-> cleavage (^-■Ccflsps-Cpray) 



br31 PARP -casp3 -> cPARP (M-M) 
br32 PARP -casp7 -> cPARP (M-M) 



brl 1* PARP -casp3 -> cPARP 
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Table 1 Summary of reactions from the original model by Bentele et al. and the reduced model (Continued) 
Inhibition of caspases-3, -7 and -9 



br33 casp3 + IAP -> casp3:IAP (MA) 


10° 


br12* casp3 + IAP -> inhibition (MA) 


br34 casp7 + IAP -> casp7:IAP (MA) 


10- 1 




br35 casp9 + IAP -> casp9:IAP (MA) 


10- 3 




br36 casp3:IAP -> casp3 + IAP (MA) 


10- 2 




br37 casp7:IAP casp7 + IAP (MA) 


10" 4 /10- 3 




br38 casp9:IAP -> casp9 + IAP (MA) 


10" 5 /10" 4 




Production of the virtual variable x aa 






br39 — > x aa (1 - x aa /Km 367act + (1 - x aa ) • (/c36act • C casp3 + k 36act ■ C casp6 + k 7act ■ 

Ccaspi)) 


10° 


brl3 * x oq fa^^ 7att - 0 - x aa yc casp3 ) 



Prototypes of parameter notations in the original Bentele's model are listed in Additional file 1: Table 1S. Abbreviations used in the table represent MA - mass action 
kinetics and M-M - Michael is-Menten kinetics. The column with rates includes rate orders for the fast/reduced activation scenarios. 



c, 



casp7' 



< a-C, 



casp3 1 



(12) 



where a = 0.18 for both activation scenarios. In addition, 
we considered inequalities holding for the parameters in 
the formula (10): 

^7 act ' Ccaspl ^ ^36act ' ^casp3 ^ ^36act ' C C asp6 1 
Km367act»l - Xaa, 

and simplified it: 

{head + 0.187c 7 ^ x 



d&aa 

dt 



Km 3 



( 1 %aa ) ' C caS p3 kdd 'X a> 



(13) 



After that we decomposed the model into modules 
(Figure 3B) corresponding to three biological steps: acti- 
vation of caspase-8, cytochrome C-induced activation of 
caspase-9 and PARP cleavage down- regulated by these 
caspases. 

Below we provide detailed description of the reduction 
process of all modules. The process is based on the flow 
chart represented in the Figure 2. If it is enough for a 
reader to have a general idea of model reduction without 
going into full details, we suggest to go directly to the 
last two paragraphs of this section, where you find the 
summary of the model reduction procedure. 

Caspase-8 activation module includes 20 species and 
14 reactions (besides the reactions of degradation), 
which could be divided into three groups: cleavage of 
procaspase-8 at the DISC complex, its activation trig- 
gered by caspase-3, and inhibition of complexes 
containing DISC by cFLIPL and cFLIPS (Table 1, 
Figure 3 A, reactions brl-brl4). 

Firstly, we eliminated slow reactions brl2 and brl3 
according to the method MR1. Next, we applied the 
quasi-steady-state analysis (MR2) to the module and re- 
moved quasi-stationary intermediates CD95R:CD95L, 



DISC:pro8, DISC:pro8 2 , DISC:p43/p41 and DISC: 
cFLIPL. Thus, in particular, we got two main reactions 
instead of brl-br6: fast formation of the DISC complex 
(Table 1, Figure 3B, brl*) and slower activation of 
caspase-8 in this complex (br2*). 

Further, we noticed that the consumptions of cFLIPL 
and cFLIPS satisfied the same kinetic laws. Therefore, 
these species could be lumped (MR3) resulting in the re- 
action br4*: 

cFLIP + 2- DISC + pro8 -> cFLIP:DISC:pro8, 

where cFLIP indicated two isoforms cFLIPL and cFLIPS 
so that C cFLI p = C cFLI pl - Ccflips- 

We approximated the concentration of caspase-6 by 
the linear function b • C casp3 , b = 0.145 (MR4), and 
merged reactions br7 and br8 of caspases-8 activation 
induced by caspases-3 and mediated by caspases-6 into 
one reaction br3*: 

pro8 -casp3 — > casp8. 

Since the reaction br8 followed the Michaelis-Menten 
kinetics with the constants k 68 and Km 68 , then the reac- 
tion rate of br3* was provided by the kinetic law 

Vbr3* — I<38 'Cp W 8'C C asp3/ (^^38 + C pro $) , 

where k 38 = 0.145 • k 68 and Km 38 = Km 68 . Analysis of 
this kinetic law for low level of CD95L ensured that 
C P ro8 » Km 3 s* In addition, in the case of the fast activa- 
tion scenario the value of v& r3 * was much lower than the 
rate of caspase-8 activation mediated by CD95. Thus, we 
established v br3 * = k 38 >C casp3 without significant changes 
in the results of the model simulation (MR1, MR5). 

Complementing the above reduction steps, we re- 
moved the elements that were unnecessary to fit the ex- 
perimental dynamics provided by Bentele et al These 
elements were degradation reactions of CD95L, CD95R 
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activation 
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activation 
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NF-kB activation 

p43-FLIP C 



casps ► 
Caspase-8 activation 

3p43-FLIP casp3<J<- 



->>casp8 

Caspase-3 activation 
— ^casp3 



Figure 4 The original model by Neumann et al. and results of the model reduction. A. The original model decomposed into modules 
according to three steps of the CD95 signaling pathways: activation of caspase-8, pro-apoptotic pathway resulting in caspase-3 activation and 
anti-apoptotic pathway regulated by NF-kB. The species retained during the model reduction are colored; reactions are represented by solid lines. 
B. The modular view of the model. Activation of caspase-8 and p43-FLIP (product of cFLIPL cleavage) occurs at the DISC complex and triggers 
simultaneous processes of cell death and survival. 
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and FLIP, and the blocked DISC complex (inactivated 
DISC). As a result, we introduced the conservation law: 

CcD95L~CcD95R = Lq-Rq = COHSt 1 

where L 0 and R 0 denote initial concentrations of the ligand 
and receptor, respectively. Solving the differential equation 



CD95R 



dt 



~ I<LR • CcD95R • CcD95L , 



we obtained the analytical function for the receptor con- 
centration: 



CcD95R 



Ro-Lo 



1 - Lq/R 0 > exp(-k L R>(Ro -L 0 )-t) 



(14) 



The module of caspase-9 activation consisted of 19 
species and 14 reactions (besides the reactions of deg- 
radation) including activation of caspase-9 triggered by 
cytochrome C, cleavage/activation of Bid by caspases-2 
and -8, IAP inhibition of caspase-9 activity and binding 
of inhibitors by Smac (brl5-br26, br35, br38). Based on 
the method MR1, we eliminated the last two reversible 
reactions, as well as the slow reaction br20 of the 
apoptosome complex dissociation. Further, since this 
complex was a quasi-stationary intermediate and the 
rate of brl5 was linearly dependent on a difference of 
brl6 and brl7 rates, we reduced the chain of 
procaspase-9 activation (brl5-br20), using the methods 
MR2 and MR4, to the form: 

Cyt C stored^Cyt C : Apaf - 1, (15) 
pro9-Cyt C : Apaf — 1— >casp9 

Accordingly to Bentele et al, we have the following 
rule of cytochrome C release from mitochondria: 



Cyt C stored 



^Cyt C stored 



(0) •) C release{t) , 



dC { 



Cyt C released 

dt 



dC( 



Cyt C stored 

dt 



Here f released) satisfies the formula (11). Taking into 
account (15) and ignoring degradation of the complex 
Cyt C:Apaf-l, we got 



dC ( 



Cyt C:Apaf-l 

dt 



dC ( 



Cyt C stored 

dt 



where c = 0.59 was the coefficient of linearity. 

Then, we noticed that the experimental measurements 
of procaspase-9 concentration were presented by Bentele 
and his colleagues only for 200 ng/ml of anti-CD95. In 



this case, degradation of procaspase-9 was insignificant. 
Thus, neglecting it and considering the kinetic law 
kApop'Cpro9'CcytC:A P af-i of the second reaction in (15), we 
obtained the differential equation of procaspase-9 dy- 



namics: 



dC, 



pro9 



dt 



-0.59.(1 ~f released) 'Q 



'^Cyt C storedif^) '^Apop 'Cpro9 



Solving it, we found: 



Cpro9 — Cprog^O)- 

f release (t) ' ( €X-\>{J<contr 't) 

exp(k contr >t) 



039 C Cyt c store d W ^Apop 



(16) 



Finally, analyzing reactions of Bid activation (br21- 
br23), we removed the slower reaction mediated by 
caspase-2. The similar reaction involving caspase-8 
followed the Michaelis-Menten kinetics with the constant 
Kwi28Bid » Csid* Therefore, we redefined the kinetics of 
this reaction based on the law of mass action (MR5): 

VbrT = k8Bid/Km8Bid-C C asp8'CBid- 

The module ofPARP inactivation contained 13 species 
(including the virtual species x aa ) and 11 reactions (ex- 
cepting reactions of degradation). Six reactions (activa- 
tion of caspases-3 and -7 by caspases-8 and -9, and 
cleavage/ inactivation of PARP) were based on the 
Michaelis-Menten kinetics (br27-br32), four reactions 
reproduced the reversible inhibition of caspases-3 and -7 
based on the law of mass action (br33, br34, br36, br37), 
and one corresponded to the production of x aa (br38), 
simplification of which was discussed above. 

We deleted slow reactions of casp3:IAP and casp7:IAP 
complexes dissociation and ignored slow degradation of 
procaspase-7 and IAP. Then analyzing the cleavage of 
procaspases-3 and -7, we eliminated reactions triggered 
by caspase-9, which were slower in comparison to the 
similar reactions induced by caspase-8. For the same 
reason we ignored the reaction of PARP inactivation by 
caspase-3 (MR1). 

Considering the remaining reactions with the 
Michaelis-Menten kinetics (br29, br32), we used the fol- 
lowing inequalities Km 78 » C pro7 and Km 376act » C PARB 
which, in combination with MR5 and (12), led to the fol- 
lowing kinetic laws: 



Vbr29 



K78 



' C caS p8 - Cpw7 1 



Km 78 

Vbr32 — t-casp3 '^PARP, 

Km 367act 

where k 3act = 0.18 • k 7act 

Thus, the CD95-signaling model consisting of 43 spe- 
cies (including one virtual species), 80 reactions and 45 
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Figure 5 Results of the composite model approximation to the experimental data by Bentele et al. The experimental data (dots) were 
obtained by Bentele et al. [17] using Western blot analysis and expressed as relative values. Once simulated values (curves in the figure) were 
found for the original Bentele's model (solid red), the reduced model (dashed black) and the composite model (solid blue), experimental points 
were recalculated taking into account these values. Wherever blue or black dots are missing, they coincide with the red dots. Note that 
considering the concentration of procaspase-8 in the original model, we observed experimental dynamics for the separate species only, but not 
for the total amount as expected from experiments. Thus, validating parameters of the composite model, we evaluated the total amount 
of procaspase-8. 



kinetic parameters was approximated by a model with 
18 species, 26 reactions and 25 parameters, except the 
constant concentration of Cytochrome C (Figure 3C). 
Figure 5 shows the comparison of simulated concentra- 
tions of these models with experimental dynamics 
obtained by Bentele et al. for the 11 species mentioned 
above. Since this dynamics was expressed in arbitrary 



units, we uniquely translated it into precise values for all 
proteins with non-zero initial values. However, for such 
molecules as caspase-8, p43/p41, tBid and cPARP, pre- 
cise values directly depended on the concentration levels 
at the end or in the middle of time series. Accordingly, 
we computed these levels during the model reduction 
and recalculated experimental dots if necessary. 
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Table 2 The role of species and reactions in the apoptotic process described by the reduced model 



Species Reactions The role in the apoptosis process 



CD95L 


br1 * 


Triggering the apoptosis process. 


DISC 


br2* 


Activation of caspase-8 in the DISC complex. In the case of a high ligand 
concentration, DISC is a slow species, so reactions br1* and br2* cannot be combined. 


caspase-3 


br3* 


Activation of caspase-8 via the negative-feedback loop. When the concentration 

of CD95L is reduced, the rate of br3* is an order of magnitude higher than the rate 

of br2*. Thus, activation of caspase-8 is mainly caused by caspase-3. 




br8* 


Cleavage/activation of procaspase-2 is confirmed by the experimental data. 




br9* 


Caspase-3 activation plays a crucial role in the cleavage of procaspases-2, -8 and PARP. 




brll* 


PARP inactivation is confirmed by the experimental measurements of PARP and cPARP 
concentrations. 


cFLIP 


br4* 


Inhibition of the DISC complex. cFLIP blocks the activity of DISC preventing a significant 
increase of caspase-8 concentration. This conclusion is consistent with the observations 
by Bentele et al. asserting that down-regulation of cFLIP resulted in cell death occurring 
already upon low concentration of anti-CD95 (1 ng/ml). 


DISC:cFLIP:pro8 


br5* 


Production of blocked p43/p41, whose concentration is detected by the experimental data. 




br6* ; br7* ; br10* 


Activation of Bid and cleavage of procaspases-7 and -9. These reactions determine 
dynamics of the mentioned species in agreement with the experimental data. 


IAP 


br12* 


Inhibition of caspase-3. In the case of the reduced activation scenario, IAP prevents caspase-3 
from reaching a significant level and, therefore, blocks significant increase of caspase-8 due to br3*. 


x aa 


br13* 


The virtual variable x aa specifies the process of degradation. 



"Experimental data" in this table refers to the data obtained by Bentele et al. [17]. 



Species without experimental evidence (CD95L, DISC, 
cFLIP, DISC:pro8:cFLIP, IAP, caspase-3 and virtual spe- 
cies x aa ) and reactions brl*-brl3* directly affect the sim- 
ulated dynamics of the experimentally measured 
concentrations (Table 2). Thus, further reduction of the 
model by removing these elements is impossible. Re- 
garding the degradation process, we agree with Bentele 
et al that experimental data cannot be matched if this 
process is fully ignored. Actually, if we remove any of 
the retained reactions of degradation, then excepted 
concentration dynamics will be impaired. Thus, we 
found the minimal complexity approximation of the ori- 
ginal model. 

Reduction of the CD95-mediated and NF-kB signaling 
model 

As was mentioned above, Neumann et al simplified 
their model in order to reduce the large number of free 
parameters. The authors noted that further simplifica- 
tion of the model was not possible because the model fit 
significantly decreased. However, we removed slow reac- 
tions (nr4, nrll-nrl3, Figure 4A) using the method 
MR1, and removed the reaction of NF-kB degradation 
(nr23), which was not significant for reproducing the ex- 
perimental data of the original model. Note that we 
retained the slow reactions nr6 and nr7, the former of 
which regulates apoptosis inhibition in the case of high 
levels of cFLIPL and procaspase-8 in accordance with 
the precursor model prediction (see the analysis of the 
predictions below), and the latter of which will be 



required to combine this model with the Bentele s model 
hereafter. 

Based on the method MR6, we also took into account 
the inequality C IKK » C p ^-flip and simplified the kinetic 
law 

V wr i9 = k p 43-FLIP IKK -Cjkk 'Cpb3-FLIP 

of the reaction nrl9 to the form 

Vnr\9* = kp43-FLIP IKK 'CiKI<(0)'C p 43-FLIP- 

Therefore, we reduced the number of model species 
from 23 to 20, the number of reactions from 23 to 18 
and the number of kinetic parameters from 17 to 15 (be- 
sides the constant concentration of IKK). These modifi- 
cations did not change the fit to the experimental data 
provided by Neumann et al (Figure 6). 

Model composition 

Analysis of the reduced models based on the Akaike cri- 
terion (7) confirmed that they had lower complexity 
than the original models (Table 3). The difference be- 
tween the mean AIC coefficients (8) of the models de- 
creased by 60%, relative to the initial value. In addition, 
the reduced Neumanns model better approximated ex- 
perimental data. Therefore, we used it as the basis for 
the model composition. 

We took into account that the considered experimen- 
tal data were obtained with SKW 6.4 [18] and HeLa [19] 
cells, which were shown to behave as type I [21] and 
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Figure 6 Comparison of the composite model simulation results with the experimental data by Neumann et al. The simulated 
concentrations of the original and reduced models (solid lines), as well as the simulated values of the composite model (dashed lines) were 
obtained by Neumann et al. [18] for three different concentrations of anti-CD95: 1500 ng/ml (black), 500 ng/ml (blue), and 250 ng/ml (red). We 
did not recalculate the relative experimental data (dots) as in the case of the Bentele's model due to a slight difference between the results of 
the composite and precursor models. An exception was made for values of procaspase-3 and caspase-3, which were normalized for the 
composite model by the term 0.1 for convenience of the visual plots comparison. 



type II [28] cells, respectively. In this regard, we assumed 
that: 

(A) initial species concentrations could vary for 
different cell lines; 

(B) kinetic parameters of all reactions (besides the 
degradation rate modeled as the function Iq • o? aa ) 
have the same values for both cell lines, whereas 
the value of Iq could be regulated by various 
entities and, moreover, is dependent on the initial 
concentration of anti-CD95 [18]; 

(C) type I and type II cells conform to different 
reaction chains of caspases-3 activation [21]. 

We constructed the composite model of CD95 and 
NF-kB signaling in the following way. Firstly, according 



to the assumption (C), we replaced reaction nrl5 in the 
Neumann s model by a chain of three reactions: 

Bid-casp8— >tBid, pro9-tBid— >casp9, pro3-casp9— >casp3 

(17) 

The first reaction of this chain had the same kinetic 
law as br7*. The other reactions represented br6* and 
nrl5 with kinetics modified using the law of mass action. 
Such modification of br6* was necessary to get the slow 
increase of caspase-3 concentration experimentally de- 
scribed by Neumann et al [19]. Modification of nrl5 
consisted in substitution of C casp9 for C casp8 . 

Secondly, we supplemented the Neumanns model 
with reaction of caspase-3 inhibition by IAP and 



Table 3 Comparison of the investigated models according to Akaike's information criterion (AIC) 





Bentele et al. 


Neumann et al. 


Composite 




Original model Reduced model 


Original model Reduced model 


model 


AIC 


182.46 138.83 


201.58 198.17 


318.10 


Mean AIC 


1.50 1.14 


0.96 0.94 


0.96 
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Table 4 Modifications of the reduced Neumann's model to reproduce the experimental data by Bentele et al 



Step Changed parameters Initial values, 
and concentrations modifications 



Reason 



Ranges and initial 
values of fitting 



1 CD95R:FADD 

2 procaspase-8 

3 FLIPS 



4 klO 

5 procaspase-3 

6 Bid, procaspase-9 



91.266 nM, increase Cose: onti-CD95 = 5 [ig/ml Procaspase-8 is cleaved in about 30 

minutes. This is possible only if DISC concentration reaches a level 
sufficient to accelerate the reaction nr2. 

64.477 nM, increase Case: anti-CD95 = 5 [ig/ml The concentration of caspase-8 should 

reach its maxim value in about 20 minutes. 

5.084 nM, increase Cose: onti-CD95 = 200 ng/ml As mentioned in Table 2, activation of 

caspase-8 is mainly caused by caspase-3 and, therefore, is delayed 
approximately for 30 minutes. Since reactions nr2 and nr5 have the 
same rates, we cannot observe the delay, but when we reduce the 
rate of nr5, we can. The latter is achieved by upregulation of the 
DISC:pro8 complex inhibition. 



0.121 nlVTmin , decrease Case: anti-CD95 
by an order of magnitude concentration. 



1.443 nM, refitting 

5.003 nM, 2.909 nM, 
replacement by 
Bentele's values 



[10 , 101,442.821, 
Bentele et al. 



Bentele et al. 

[10 1 , 10 2 ], 65.021, 
Bentele et al. 



-- 200 ng/ml Preventing a rapid growth of caspase-8 0.012 



Case: anti-CD95 = 5 \ig/ml Improving approximation of the 
experimental data by Bentele et al. 

Case: anti-CD95 = 200 ng/ml Reproducing the experimental dynamics 
of procaspase-9. 



[10°, 10 1 ], 1.443, 
Neumann et al. 

231.760 (Bid), 
245.101 (pro-9), 
Bentele et al. 



evaluated parameters in this reaction together with 
parameters in (17) using the corresponding data by 
Neumann et al and optimization tools of BioUML [29] . 

Thirdly, we analyzed the changes of the initial concen- 
trations and parameters of the derived model required 
to reproduce the Bentele s experimental data fixing levels 



of procaspases-3, -8, cleaved product p43/p41, and 
caspase-8 (Table 4). All the changes, except step 4, 
agreed with assumptions (A) and (B). To reproduce the 
species dynamics, we had to increase the initial concen- 
tration of procaspase-3 in the case of HeLa cells by an 
order of magnitude. 
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Figure 7 Creation of the composite model. The modules of caspase-8 and NF-kB activation were taken from the reduced Neumann's model 
and combined with the modules of caspase-9 activation and PARP inactivation isolated from the reduced Bentele's model. In addition, two last 
modules were supplemented by reactions from the modified module of caspases-3 activation. 
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Finally, we combined the modules of caspase-8 and 
NF-kB activation with the modules of caspase-9 acti- 
vation and PARP inactivation (Figure 7). For this pur- 
pose, we modified the last two modules based on (17) 
and refitted their parameters as well as the concentra- 
tions mentioned in Table 4, using experimental data 
by Bentele et al. For better fit we also supplemented 



the model by degradation reactions of some species 
(procaspase-9, caspase-9, DISC:pro8 and DISC:pro8: 
FLIPS). 

The resulting composite model (Figure 8) consists of 
30 species (including one virtual species), 38 reactions 
(including 15 reactions of the species degradation) and 30 
kinetic parameters, except for the constant concentration 




ApoptOtiC activity 
Caspasc-8 activation 

P43-FUP casp8 casp3<j<- 

-» 
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NF-kB activation 



->>Apoptctic activity 
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Figure 8 The composite model of the CD95- and NF-KB-signaling. A. The model integrating pro- and anti-apoptotic machinery described by 
the models by Neumann et al. and Bentele et al. Reactions with the "nr" titles were taken from the first model, while the "br" titles indicate 
reactions from the second one. The "*" symbol marks reactions defined in the reduced models. The subscript "m" means that reaction was 
modified for the models composition. B. The modular view of the model. Activation of caspase-8 triggered by CD95 leads to NF-kB activation 
and cell survival, on the one hand, and PARP cleavage resulted in the cell death, on the other hand. 
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of IKK (Additional file 1: Table 2S, Additional file 1: Ta- 
ble3S, Additional file 1: Table4S). The model showed a 
good fit between the simulation results and the experi- 
mental data (Figures 5 and 6). 

Calculation of the mean AIC coefficient for this model 
revealed that it has the same level of complexity as the 
reduced models (Table 3). We also noted that reduction 
of both models did not significantly alter the steady-state 
concentrations (Additional file 1: Table 5S, Additional 
file 1: Table6S). However, the composite model has dif- 
ferent steady-states, which are more stable according to 
the mean sensitivities (9) in the cases of SKW 6.4 cells 
and HeLa cells stimulated with 250 ng/ml of anti-CD95 
(Additional file 1: Table 7S). 

The reduced models are equivalent to the precursor 
models with respect to their biological predictions 

Analyzing the models constructed by Neumann et al 
and Bentele et al, we divided predictions made by them 
into two groups: qualitative and quantitative. The first 
group concerned the model network and described me- 
chanisms of the protein-protein interactions. It is this 
group to which we assigned the experimentally validated 
Neumann's predictions that processed caspases are not 
required for NF-kB activation, as well as that pro- 
apoptotic and NF-kB pathways diverge already at DISC. 
The corresponding network of reactions was preserved 
in the reduced model. Thus, the prediction remained 
valid. 

The second group of predictions characterized behav- 
ior of the models after changing concentrations of some 
species, such as CD95L, CD95R, procaspase-8, and in- 
hibitors cFLIPL, cFLIPS, and IAP. Analyzing the predic- 
tions of the models by Neumann et al and Bentele et al 
listed in the Tables 5 and 6 respectively, we concluded 
that the reduction of the first of them did not alter the 
species dynamics (Additional file 1: Table 8S). For the 
second reduced model, we detected some minor changes 
in the predictions. However, in general, the model be- 
havior remained in a good agreement with the Bentele's 
experiments (Additional file 1: Table 9S). 

Model composition modifies some properties of the 
precursor models and leads to new predictions 

Considering the predictive ability of the composite 
model, we found that in the case of HeLa cells it was 
completely preserved (Table 5). 

In the case of SKW 6.4 cell line, we estimated the 
value of degradation rate parameter kd based on the ex- 
perimental observations by Bentele et al for 1-10 ng/ml 
of anti-CD95 (Table 6). When performing simulation 
experiments of the model for this cell line, we used 
three different values of kd depending on whether the 
initial concentration of CD95L was within the range 



1-100 ng/ml, within the range 100-1000 ng/ml or 
greater than 1000 ng/ml (Additional file 1: Table 3S). 

Analyzing the model predictions, we detected differ- 
ences in the behavior of the Bentele s and the composite 
models (Table 6). The first of them describes a threshold 
mechanism for CD95-induced apoptosis implying that 
this process is completely stopped when CD95L concen- 
tration is below the critical value. However, the model 
predicts a sharp increase of cPARP immediately when 
the level of CD95L exceeds the threshold (Table 6, Ne 1). 
In contrast, the composite model shows a gradual in- 
crease of cPARP with increase of the ligand level. To de- 
termine the apoptotic threshold in this case, we found 
concentration of CD95L for which cPARP ratio is equal 
to 10% of the initial PARP amount [30]. As in the ori- 
ginal model, this concentration was in the range of 1-10 
ng/ml. Additionally, the models revealed sensitivity of 
the threshold to the concentration of cFLIP (Table 6, 
2). Both of them predicted the cell death phenotype 
upon stimulation by 1 ng/ml anti-CD95 and the level of 
cFLIP decreased by -50%. 

The next prediction, which we observed for SKW 6.4 
cells, considers a delay of caspase-8 activation caused by 
low concentrations of anti-CD95 (Table 6, 3). The 
maximal delay predicted by the Bentele's model was 
about 2500 min, whereas the composite model detected 
a much lower value (700-800 min). Analysis of depend- 
ence of caspase-3 concentration on initial values of IAP 
and CD95L (Table 6, 4) also demonstrated differences 
in the models' predictions. Thus, low level of IAP (less 
than 1 nM) in the original model results in complete cell 
death. However, under the same condition the compos- 
ite model shows a smooth increase of caspase-3 levels 
with CD95L changing from 0.3 nM to 1000 nM. If the 
amount of the ligand is less than 0.3 nM, IAP blocks 
apoptosis completely. In addition, the Bentele's model 
predicts that high concentration of IAP prevents a sig- 
nificant increase of caspase-3 even for high levels of the 
ligand, whereas in the composite model high concentra- 
tion of CD95L leads to cell death. 

Limitations of the composite model 

We analyzed a set of parametric constraints which 
allowed us to reduce the precursor models (Additional 
file 1: Table 10S). Following them, we can recover a de- 
tailed description of the investigated process. Connec- 
tion of the composite model with the original model by 
Neumann et al was preserved, whereas with the original 
model by Bentele et al it was partially lost since the 
chains of species interactions were changed for the sake 
of combining the models together. In particular, the 
pathway of caspase-8 activation in DISC was derived 
from the Neumann's model as well as the reaction 
chain of procaspase-9 cleavage was modified. However, 
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Table 5 Analysis of predictions regarding apoptosis in HeLa cells as formulated by Neumann et al 
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The concentration of anti-CD95 required for the apoptosis induction (the apoptotic 
threshold), is within the range of 30-100 ng/ml. This range remains the same for CD95 
decreased by about 12-fold. The simulation time, which we used to reproduce this 
prediction, was 60 hours. 



0 50 100 150 

anti-CD95 cone, ng/ml 




500 1000 

Time, min 



The decreased receptor number results in impairment of both CD95- and NF-KB-signaling 
pathways. To test this prediction, Neumann et al. considered levels of caspases-8 cleavage 
and kB-a degradation for the original (solid lines) amount of CD95 and the amount 
decreased by about 12-fold (dashed lines). The concentration of anti-CD95 was 500 ng/ml. 
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Along with increasing the concentration of anti-CD95 from 500 ng/ml to 1500 ng/ml, p43/ 
p41 peaks earlier, while there is almost no difference for p43-FLIP. 
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Increased concentrations of cFLIPS inhibit both apoptotic and NF-kB pathways, although 
p43-FLIP generation is inhibited at a lower threshold than p43/p41 generation. 
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Increasing the concentration of cFLIPL leads to a steep increase in p43-FLIP generation until 
it reaches a maximum, after which the curve drops. Lowered levels of cFLIPL result in very 
little p43-FLIP but almost unchanged levels of p43/p41. 

At very high concentrations of cFLIPL no p43-FLIP is generated. This drop-off was not 
observed experimentally by the authors. 
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Table 5 Analysis of predictions regarding apoptosis in HeLa cells as formulated by Neumann et al (Continued) 




Only an intermediate level of cFLIPL promotes NF-kB activation. Decreased levels of 
procaspase-8 lead to a significantly lower amount of p43-FLIP and, subsequently, NF-KB.The 
figures show of logarithmic dependence of the maximal NF-kB concentration on the initial 
values of procaspase-8 and cFLIPL 
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Procaspase-8 init.conc, nM 




High cFLIPL or low procaspase-8 concentrations cause suppression of apoptosis. The figures 
show the same dependence as considered in the previous prediction, but with caspases-3 
instead of NF-kB. 
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All the predictions marked with an asterisk were experimentally tested by Neumann et al. and confirmed, unless otherwise noted. The simulation time in 
predictions 3-7 was 360 min. The concentration of anti-CD95 considered by the authors in predictions 4-7 was 1000 ng/ml. 



following the constraints of Additional file 1: Table 10S, 
we can complete the composite model with such reactions 
as Bid truncation by caspase-2 and PARP cleavage trig- 
gered by caspase-7. We can also recover the details of 
caspase-9 activation mediated by cytochrome C and Smac 
release from mitochondria but with altered kinetics of re- 
lease. For now, the question of admissibility of such modi- 
fications as well as extension of the composite model 
based on other apoptosis models remains open and is a 
challenge for further research. 

Discussion 

In this paper, we considered two approaches to develop- 
ment of mathematical models of cell fate decisions. The 
first concerns the methodology of model reduction and 
involves approximation of one model by another one of 
lower dimension without affecting dynamics of experi- 
mentally measured species. The second implies compos- 
ition of the models and aims at reproducing experimental 
dynamics of all precursor models. 

There are many advantages of working with low- 
dimensional models [8]. In particular, the researcher has 
a clear vision of the most important biochemical reac- 
tions taking place in the modeled system as well as bet- 
ter understanding of them. Low-dimensional models are 
easier to analyze and faster to simulate. This helps to 
save time and enhances productivity. The main limi- 
tation of model reduction consists in loss of biological 
information. However, it should be noted that even if 



some information (for example, about very slow bio- 
chemical reactions) may be lost, it can result in a more 
clear understanding of the most important interactions 
and allows focusing on the decisive processes in the 
model, predictive ability of which is reasonably pre- 
served. Hence, this limitation can be transformed into 
an advantage. 

Model composition aiming at getting a single model 
from several ones is useful because in such a case a 
computational biologist is able to investigate the com- 
posite model behavior under different conditions that 
cannot be performed in the precursor models separately. 
For example, our model allows studying the role of p43- 
FLIP or IAP in the type I SKW 6.4 or type II HeLa cells, 
respectively, that might become a task for the future 
work. In other words, we constructed the model that de- 
scribes pro- and anti-apoptotic signal transduction in 
different cell types with reasonable accuracy instead of a 
couple of different models. Whenever necessary, some 
reaction chains and parameters can be switched off giv- 
ing opportunity to simulate a certain type of cells. In 
addition, the composite model covers experimental data 
obtained from all precursor models, each of which sep- 
arately satisfies its own data only. 

Model composition sometimes causes modification of 
some properties of the initial models, resulting in new 
testable predictions. In our case, such predictions were 
related to SKW 6.4 cells, and some simulation results 
were different from the corresponding results of the 



Kutumova et a I. BMC Systems Biology 201 3, 7:1 3 
http://www.biomedcentral.eom/1 752-0509/7/1 3 



Page 19 of 21 



Table 6 Predictions of the models for SKW 6.4 cells 
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Experimental observations by Bentele et al. and predictions of the models 
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Experiments by Bentele et al.: 

- for 1 ng/ml of anti-CD95, PARP cleavage was not observed; 

- the measured death rate for 10 ng/ml of anti-CD95 was 20-30%. 
Original model (red): 

- the apoptotic threshold is 1.9 ng/ml; 

- cPARP concentration rises dramatically within an extremely narrow interval of anti-CD95 
levels overcoming the apoptotic threshold. 

Composite model (blue): 

- the apoptotic threshold is 3.5 ng/ml; 

- cPARP concentration rises in a smooth manner along with the increase of anti-CD95 level. 

Experiments by Bentele et al.: down-regulation of cFLIP in SKW 6.4 cells by addition of 
cyclohexamide resulted in cell death (40% for 1 day) already upon 1 ng/ml of anti-CD95. 
The level of cFLIP was decreased to 70%. 

Original (red) and composite (blue) models: 

- the apoptotic threshold is highly sensitive to the concentration of cFLIP; 

- decreasing the initial concentration of cFLIP by more than 51% and 49% for the original and 
composite models, respectively, leads to cell death upon stimulation by 1 ng/ml of anti-CD95. 



2 500 
| 2 000 
^ 1 500 
^ 1 000 



TO 



500 
0 




0 5 10 15 20 25 30 

anthCD95 cone, ng/ml 



Experiments by Bentele et al.: in the 10 ng/ml activation scenario, a significant increase of 
caspase-8 was observed after more than 4 hours. 

Original (red) and composite (blue) models: anti-CD95 concentrations which are slightly 
above the apoptotic threshold result in caspase-8 activation after a delay of many hours. 

The figure shows peak times of caspase-8 concentration exceeding 0.1% of the initial 
procaspase-8 level. 
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Original model: 

- low concentrations of IAP (less than 1 nM) result in complete cell death; 

- high concentrations of IAP prevent a significant increase of caspase-3 even 
for high concentrations of the ligand. 

Composite model: 

- low concentrations of IAP (less than 1 nM) block apoptosis for CD95L less than 0.3 nM; 

- high concentrations of CD95L lead to cell death. 

The figures show logarithmic dependence of the maximal caspases-3 concentration 
on initial values of IAP and CD95L. 
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The simulation time was 2880 min (2 days) in all predictions. The apoptotic threshold in the prediction 1 is the concentration of anti-CD95 after which cPARP 
amount exceeds 10% of the initial PARP level. 
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original Benteles model Nonetheless, the composite 
model behavior remained in good agreement with ex- 
perimental data used in modeling by Bentele et al. 
Thus, the question of what model (original or com- 
posite) is more realistic requires further experimental 
investigation. 

One of the questions that can be asked concerns the 
need of reducing models in order to combine them. An 
attempt to construct the composite modular model of 
apoptosis was made in [26]. This attempt revealed mul- 
tiple conceptual and methodological difficulties. The 
main obvious obstacles included: 

- choice of elements in the overlapping parts of 
different models; 

- incomparable sets of parameters; 

- the lack of experimental data, as well as inability to 
use the data obtained for different cell lines or in 
different ways (e.g. single-cell or cell culture 
measurements); 

- inability to make accurate predictions. 

Model reduction allows solving some of these prob- 
lems. In particular, reducing the number of model ele- 
ments, we reduce the overlapping parts as well. This 
may be essential for direct combining of models. In 
addition, the complexity of the reduced model become 
comparable with the complexity of available experimen- 
tal data. Therefore, the risk of model overfitting is 
decreased. 

In the case when two models are directly related (for 
example, one model was emerged from the other), their 
composition may be significantly easier in comparison 
to composition of quite different models. In our work, 
the models by Bentele et al. and Neumann et al. are not 
directly related, as it could be expected as these models 
were constructed in the same research group. For ex- 
ample, they use different reaction chains of caspase-8 ac- 
tivation and different values of kinetic parameters in the 
overlapping reactions. 

Another quite useful principle that we used in model- 
ing was modular structure of the developed model. This 
principle provides flexibility for future extensions. Thus, 
we are planning to extend the composite model and 
supply it with modules and data from studies of TRAIL 
[31-33] and TNF-a [34,35] signaling, apoptosome- 
dependent caspase activation [36] and p53 oscillation 
system [37]. It is noteworthy that the model by 
Laussmann et al. [33], describing TRAIL-induced ac- 
tivation of caspase-8, emerged from the original study 
by Bentele et al. [18]. This fact may help us in our 
future work. 

Characterising the models used in our work, we can 
say that we saved our time working with the model by 



Neumann et al. [19] derived from the BioModels data- 
base [27], and spent a lot of time to reproduce the 
model by Bentele et al. [18] that had not existed in a 
ready-to-use format. Thus, we would like to emphasize 
that for the common benefit, the systematic application 
of biological standards in modeling (e.g. SBML [38] and 
SBGN [39]) and depositing the working models in public 
databases (e.g. BioModels [27]) would significantly facili- 
tate analysis of existing biomathematical models and 
using them as the base for development of new compos- 
ite systems. 

Finally, we believe that the composite model itself may 
be useful for further investigation of apoptosis. 

Conclusions 

Mathematical modeling provides a powerful tool for 
studying the properties of biological processes. Methods 
of model reduction allowed us to take a first step to- 
wards validation of the modular model of apoptosis [26]. 
Using these methods, we composed two models describ- 
ing pathways of CD95- and NF-i<B-signaling in one 
without affecting the fit to the experimentally measured 
dynamics and model predictions. For the model reduc- 
tion and composition, we used the BioUML software 
that was extended by the required methods of analysis. 

The models by Bentele et al. and Neumann et al. 
reconstructed in BioUML and represented using SBML 
format and SBGN notation, as well as their reduced and 
decomposed versions are available at http://ie.biouml. 
org/bioumlweb/#de=databases/The%20composite%20mo 
del%20of%20CD95%20and%20NF-kB%20signaling. 
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